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A new measure to characterize stability of complex dynamical systems against large perturbation 
is suggested, the stability threshold (ST). It quantifies the magnitude of the weakest perturbation 
capable to disrupt the system and switch it to an undesired dynamical regime. In the phase space, the 
stability threshold corresponds to the “thinnest site” of the attraction basin and therefore indicates 
the most “dangerous” direction of perturbations. We introduce a computational algorithm for 
quantification of the stability threshold and demonstrate that the suggested approach is effective 
and provides important insights. The generality of the obtained results defines their vast potential 
for application in such fields as engineering, neuroscience, power grids, Earth science and many 
others where robustness of complex systems is studied. 



2 



Figure 1 . Stability threshold (ST) and its quantification. Attractor A, its attraction basin B and ST a. The trace of the 
algorithm converging to the point M is shown by black dots. Other LOOT points M2 and M3 are also shown. In the the 
zoomed part, safe and unsafe perturbations are shown. Dotted lines are perturbations, solid black lines are trajectories of the 
perturbed system. 


INTRODUCTION 

Complex systems science is strongly based on linear stability analysis considering small perturbations of dynamical 
systems. In a seminal paper [T] this concept was extended even to the stability of synchronization in complex 
networks leading to the efficient master stability formalism. However, for various applications often the influence of 
large perturbations is also of crucial importance. Typical examples are climatological systems, in particular ocean 
circulations. Well accepted is that the Atlantic Meridional Overturning Circulation may be sensitive to changes 
in the freshwater balance of the northern North Atlantic. When an anomalous freshwater flux is applied in the 
subpolar North Atlantic, this circulation collapses in many ocean-climate models [2]. Another example is power 
grids which are networks of connected generators and consumers of electrical power. For proper function of such 
networks synchronization between all the nodes is essential. Local failures, overloads or lines breaks may cause 
desynchronization of nodes and lead to large-scale blackouts mm- 

The study of system’s stability against large perturbations implies treating the following challenging problem: 
definition of the class of “safe”, or admissible perturbations after which the system returns back to the initial regime. 
In contrast, “unsafe” perturbations switch the system to another, often unwanted, dynamical regime. The definition 
of the class of safe perturbations of a nonlinear system is very complicated and basically different from the linear 
stability analysis. The reason is that for large perturbations linearization is inadequate and the perturbed dynamics 
is governed by nonlinear equations whose analytical study is impossible in general. Some analytical methods do exist, 
for example the method of Lyapunov functions [5]. However, this method has serious limitations since a Lyapunov 
function for a particular dynamical system is often not constructive. The “safe” perturbation class can also be 
analytically estimated for some specific systems, e.g. networks of spiking neurons 0. Nevertheless, an important task 
is to develop numerical methods of defining and describing the class of safe perturbations [7] . 

From the viewpoint of nonlinear dynamics, established dynamical regimes of the system corresponds to attractors 
in the phase space. The class of safe perturbations is equal to the attraction basin, i.e. the set of the points which 
converge to the attractor. A perturbation is safe if it brings the system to a point within the attraction basin. The 
first attempt to characterize attraction basins in complex networks was undertaken in |8j where the concept of basin 
stability was introduced. The basin stability equals 

Sb= x(x)p(x)dx, (1) 

Jq 

where Q is the set of possible perturbed states x, p{x ) with Jq p(x)dx = 1 is the density of the perturbed states, 
and x( x ) equals one if the point x converges to the attractor and zero otherwise. The value Sb £ ( 0 ; 1 ] expresses the 
likelihood that the perturbed system returns to the attractor. An important advantage of this measure is that it can 
be easily calculated by Monte-Carlo method. Namely, one should just take a large number of random points in the 
phase space and check how many of them converge to the attractor. If M of N 3> 1 points converge, Sb ~ M/N. 

Basin stability is an important characteristic extending the concept of linear stability for the case of large perturba¬ 
tions. However, many real dynamical systems, especially complex networks, possess highly-dimensional phase space 
with complicated structure. This makes it problematic to characterize an attraction basin by just a single scalar 
value. Moreover, basin stability depends on the perturbation class Q which should be chosen a priory. 
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In this paper we suggest a new measure to characterize stability against large perturbation, the stability threshold 
(ST). We were inspired by the observation that for real systems it is often important to know the maximal magnitude 
of perturbation which the system is guaranteed to withstand, like the maximal voltage jump for a stabilizer or the 
maximal bullet energy for a bulletproof vest. In the following, in Sec. 1 we introduce the ST in detail and explain 
how to calculate it. In Sec. 2 and 3 its potential is demonstrated for two paradigmatic model systems. In Sec. 4 we 
relate the two measures, the stability threshold and the basin stability. In Sec. 5 we briefly summarize and discuss 
our results. 


I. DEFINITION AND QUANTIFICATION OF STABILITY THRESHOLD 

We define ST as the minimal magnitude of a perturbation capable to disrupt the established dynamical regime, i.e. 
to push the system out of the attraction basin. In the phase space, ST is the minimal distance between the attractor 
A and the border 6B of its attraction basin, i.e. 


a = inf {dist(a, b) |a G A, b G 5B} , (2) 

where dist(-, •) is the distance between the points in the phase space. Further, we use Euclidean metric to calculate 
the distance. However, any other metrics can be used as well. 

To better understand the physical meaning of ST consider the system settled to the attractor A as depicted in Fig. 
1. Let a € A and b G SB be points corresponding to ST such that dist (a, b) = a. Consider now a perturbation of 
the system which results in a quick change of its state. Namely, let the system state change from the unperturbed 
one x u to the perturbed one x p so that the perturbation Aa’ = x p — x u has the magnitude q = |Ax| . If q < a, the 
perturbation can never kick the system out of the attraction basin (Aaq in Fig. 1). But if q > a and the system is 
near the point a just before the perturbation, it may be kicked out of the basin if the direction of the vector Aa: is 
close to the vector D = b — a (Ax 2 in Fig. 1). The above reasoning shows that besides the value of a, the direction 
of the corresponding vector D is critical. This vector corresponds to the most “dangerous” direction of perturbations 
in which the distance to the basin border is the shortest. 

Let us now come to a question of quantification of the ST. For this purpose we suggest a two-stage algorithm whose 
basic principles are described below and also illustrated in Fig. 1. 

i) On the first stage we identify some point K\ on the border of the attraction basin. For this purpose we choose 
an arbitrary point Kq in the vicinity of the attractor and start to move from the attractor until leaving the basin. 
The point K\ is found then by the bisection method. 

ii) On the second stage we move along the basin border. On each step we draw a tangential hyperplane to the 
border at the current point K n . In the hyperplane we find the point closest to the attractor A and make a step 
towards this point and so obtain the new point K n +i- Such steps bring us closer and closer to the attractor and 
finally converge to the point M on the border with the minimal distance to the attractor . 

The second stage relies on smoothness of the border for only in this case it can be approximated by a tangential 
hyperplane. Further we assume the border is smooth (defined by a function of class C 1 ) which is the case for many 
realistic systems. 

The suggested algorithm allows us to determine the local minima of the distance between the attractor and the basin 
border, which we call further “local threshold” (LOOT) points. Starting from different initial points we get different 
LOCT points Mi, M2, ...,M m (Fig. 1). Between them, the one closest to the attractor is the “global threshold point” 
corresponding to the ST: a = min(cri,..., a m ), where aj = dist ( Mj,A ). 

This brute-force search to obtain all the local minima does not seem to be a very effective strategy. However, 
effectiveness of the method is essentially improved in a parametric study, i.e. when the system properties are studied 
versus its parameters. Note that such tasks are typical since all realistic systems depend on parameters and usually 
one wants to know what happens if they are varied. Suppose that for a certain parameter value p = po we have found 
all LOCT points Mi(p 0 ),M m (p 0 ). In a robust system, the phase space structure changes continuously when p is 
changed. Thus, the coordinates of LOCT points depend continuously on p. So, when p is changed by a small value 
A p = p — p 0 one should start the algorithm from the points Mj(p 0 ). Since the actual positions of Mj(p) are close, 
the algorithm converges to them quickly. In this manner one can effectively trace the positions of LOCT points over 
the parameter value. 

One can still argue that to find all the LOCT points even for one parameter value po may be a practically impossible 
task for complex high-dimension systems. However, it is important to emphasize that this task is significantly 
simplified if the system is a network, i.e. it is composed of many low-dimensional subsystems. In this case the 
coupling strength is a natural choice for the parameter p one can trace the system over. For the case of no coupling 
(p = 0) the high-dimension phase space of the whole system is just a direct product of the low-dimension subspaces. 
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Figure 2. Stability of the pendulum (|3|. (a) Phase space for a = 0.04, P = 0.1. Green area is the attraction basin of the steady 
state O, red curve is the limit cycle L. LOOT points are depicted by red dots, traces of the algorithm by black, (b) Local 
stability thresholds cti (red), 02 (blue), and <73 (black) versus P. (c) The basin stability values Sb 1 (red), Sb2 (blue) and Sbo 
(black) versus P, mean values and variances. 


In each of these subspaces the LOOT points can be found relatively easily which allows to spot them in the full space 
as well. Further one can gradually increase the coupling and trace the points. 

Another issue is the possible emergence of new LOOT points. We have an effective method to trace once found 
points over the parameter, but how can one assure that no other points have emerged closer to the attractor? This 
problem is typical for the global extremum seeking in the case when full sampling of the space is impossible [9]. 
Although there is no way to guarantee that the found ST point is indeed the global minimum, there is a way to make 
the probability of that arbitrarily close to one. For this sake one has to check sufficiently many random initial points 
and make sure that none of them gives a better result. We revisit this issue in the end of Sec. 4 where we provide an 
estimate of the number of trials one should run to decrease the mistake probability below a given level. 


II. STABILITY THRESHOLD OF PENDULUM 

Now we show how the ST approach can be applied to study some paradigmatic dynamical systems. First we 
consider a classic pendulum under an external force P: 


— =uj, —— = — auj + P — sin 61 (3) 

at at 

Here, 6 and w are the deviation angle and the angular velocity, and a describes friction. Noteworthy models similar 
to ([3]) are often used to describe dynamics of nodes of power grids, i.e. generators or consumers [4l ITf)]. The phase 
space of the model is a cylinder S' 1 x R 1 and includes two attractors: a stable steady state 0(arcsinP, 0) and a stable 
limit cycle L (Fig. 2a). In the context of power grids, the steady state corresponds to the state when the generator 
operates in synchrony with the grid, and the limit cycle corresponds to an undesired asynchronous regime. 

Next we use the concept of ST to study the attraction basin of the steady state O. The identified LOOT points are 
depicted by red dots in Fig. 2a. The most important ones are M\ corresponding to positive perturbations and M 2 
corresponding to negative ones. Because of the complex shape of the attraction basin other LOOT points exist further 
from the attractor, e.g. M 3 . Figure 2b demonstrates the local thresholds (LOCTs) Uj = dist (0,Mj) associated with 
these points in dependence on the parameter P. One can see that for P < P* « 0.15 the point closest to the attractor 
is Mi, while for P > P* the closest point is M 2 ■ Thus, the ST equals a± for P < P* and <72 for P > P *, i.e. the most 
dangerous are positive perturbations for small P but negative perturbations for large P. 

It is interesting to compare both basin measures: stability threshold and basin stability. For this sake Sb is plotted 
versus P in Fig. 2c. We calculate it for three different classes of perturbations: positive perturbations (Sbi for 
Q 1 = [—7r; 7r] x [0; 3]), negative perturbations ( Sb 2 for Q 2 = [—7r; 7r] x [—3; 0]), and perturbations of both signs ( Sbo 
for Q = Q 1 U Q 2 )- When P increases, the basin stability for all classes of perturbations decreases as well as the ST. 
Thus, both measures indicate that the system becomes less robust. However, basin stability fails to detect which 
perturbations are more dangerous: Sb 2 is sufficiently larger than Sbi for all values of P. 
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We also checked that the efficiency of our algorithm is essentially improved by tracing LOCT points over the 
parameter. For this sake we identified the position of the point Mi for P £ [0.1; 0.4] for three different setups: 
parameter step A P = 0.02 (16 datapoints), without tracing; the same parameter step, with tracing; smaller step 
A P = 0.01 (31 datapoints), with tracing. Without tracing, the search started each time from the same point. With 
tracing, the search for the new parameter value started from the position found for the previous parameter value. The 
total computation time T c equals 9 x 10~ 3 (a.u.) for the first setup, 19 x 10 -4 for the second setup, and 24 x 10 -4 for 
the third setup. Thus, with tracing the computation time decreases approximately five times. For higher-dimensional 
systems the improvement is even much higher. Note also that T c in the third setup increases by less than 30% with 
respect to the second setup, although the number of datapoints is twice larger. The reason is that with a smaller 
parameter step the positions of the LOCT points change less and they are found faster. 


III. STABILITY THRESHOLD OF COMPLEX NETWORKS 

The second example is a network of coupled one-dimensional maps. We chose maps for two reasons: first, because 
of simpler implementation, and second, to demonstrate the generality of our approach. The network on N nodes is 
governed as follows: 


N 

Xi(t + 1) = axi(t) + bxf(t) + K^Cij ( Xj(t) - Xi(t)). (4) 

1=1 

Here, 0 < a < 1 is the system parameter, k stands for the global coupling coefficient and Cij are the elements of 
the coupling matrix. Coupling between two nodes i and j equals KCij. The network has the only attractor, the stable 
fixed point 0(0,0, ...,0). However, after a large perturbation the system trajectories may go to infinity. 

For network Q, a natural way to find LOCT points is to trace them over the coupling coefficient k. For k = 0, 
the nodes are uncoupled and each of them is governed by the map Xi(t + 1) = axi(t) + xf(t), which has a stable 
fixed point Xi = 0 with the attraction basin — 1 < Xi < 1 — a. The borders of this interval define two LOCT points 
in the network phase space: Mi + (xj = 0 (j ^ i), Xi = 1 — a) corresponds to positive perturbation of the node *, and 
Mi _ (Xj = 0 (j ^ i), Xi = —1) to negative ones. We start from these points for k = 0, then gradually increase k and 
trace their positions. We also periodically check for emergence of new LOCT points, but failed to detect any. 

We study various networks with 2 < N < 100 and different types of topology: all-to-all, random nu, small-world 
P2|, scale-free m , and cluster networks HU- In all the cases, the behavior of LOCT points is quite similar. When k 
increases, the positions of the points change, so that the coordinates Xj (j ^ i ) of Mi± are no longer zeros. However, 
for weak coupling the coordinates of LOCT points obey |aq| ^$> Xj , i.e. the corresponding perturbation mainly concerns 
the node i. For larger k the situation changes and LOCT points may have several coordinates of the same order. 
Typical LOCT points are illustrated in Fig. 3a. 

Now let us consider LOCTs cq± = dist(0,_Mj±) associated with LOCT points. A typical dependence of these 
thresholds on k is illustrated in Fig. 3b. For all j, cq + grows with k, while <Ti_ decreases. Some of the points Mj + may 
disappear at certain k as well. A detailed study shows a remarkable feature of the LOCTs eq±: they turn out to be 
strongly correlated with the values of total connections strength to the node Ki = Cij- In Fig- 3c the LOCTs 

Oi± are plotted versus Ki for various nodes, coupling coefficients, network sizes and configurations. The correlation is 
large, especially for small Ki. Notice that for Ki < k* « 0.6 positive perturbations have a lower threshold than negative 
ones, and this threshold increases with Ki. This finding leads to an easy and intuitively clear rule: the stronger the 
node is connected to the network the harder it is to tear it off. However, too strong coupling (k^ > k*) is undesirable, 
since it increases susceptibility to negative perturbations. 

The global ST of the network is defined by the lowest LOCT. Figure 3b illustrates a typical dependence of the ST 
on k. 


IV. STABILITY THRESHOLD AND BASIN STABILITY 

It is interesting to compare the two measures of stability against large perturbations, the stability threshold and 
the basin stability for the same network (Fig. 3d). As the perturbation class Q we use a hypersphere of radius q = 0.8 
with constant density p which means that we consider perturbations of amplitude q and random direction. |15| One 
may see that Sb = 1 when the ST exceeds q for k > K q ~ 0.26. This confirms that the ST indeed characterizes the 
weakest perturbation that can disrupt the network. 
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Figure 3. Stability of the network Q. (a) Coordinates of a typical LOOT point M;_|_ (i = 5) for two values of k - small (blue) 
and large (red), (b) LOCTs Ui± versus k. Blue thin curves for cr; + , green thin curves for er^-. Red thick curve is the ST a. (c) 
LOCTs <Ji± versus nodal coupling strength m for various nodes, coupling coefficients and network configurations. Blue dots 
for <Ji+, green for er;_ (d) The basin stability versus k. The inset shows the same dependency in logarithmic scale. Red curve 
for numerical results, blue thin line for the estimate 



P(B) 


Figure 4. (a) Estimate of basin stability for perturbations slightly exceeding the ST. (b) The posterior probability P(B\Mi) 
versus the prior probability P(B) for various values of a denoted on the figure. 



From Fig. 3d one may acquire the wrong impression that the basin stability reaches unity much earlier than n 
reaches n q . The reason is that Sb approaches unity very quickly when a approaches q. This can be seen in the inset 
of Fig. 3d which has a logarithmic scale. This feature seems to be typical for high-dimensional dynamical systems. 
Indeed, consider an arbitrary dynamical system in the ./V-dimensional phase space settled into the attractor A with 
the stability threshold a. Consider the perturbation class Q consisting of perturbations with the amplitude q. For 
q < cr, the set Q resides inside the attraction basin B , therefore Sb = 1. For cr = q, the set Q contacts the border of 
the basin SB. For q > a some part of the set Q gets out of the basin B and Sb becomes smaller than one (Fig. 4a). 
The probability of the perturbed state to be out of the basin is proportional to the surface area s of the protrusive 
part (gray in the figure), so 1 — Sb ~ s. To estimate the surface area, one can approximate both surfaces Q and SB by 
quadratic forms near the site of their intersection. Then, the transverse size of the protrusive part can be estimated 
as d ~ sjq — cr, and the surface area s ~ d^ -1 . This leads to the estimate 


1 - S B ~ (q- c^Y 


(5) 


The corresponding slope is given by the blue line in the inset of Fig. 4d and agrees with the numerical results. 

The scaling law ([5]) suggests that for high-dimensional systems it is very unlikely that the system will be disrupted 
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by a random perturbation of magnitude just above the ST. From the other side, a wisely designed perturbation can 
disrupt the system even being slightly above the ST. 

From computational prospective, the estimate <© shows that attempts to estimate the stability threshold from 
basin stability is inefficient for high-dimensional systems. Indeed, on the example from Fig. 4d one can see that 
it is very complicated to detect the exact point where the basin stability reaches exactly unity. The probability to 
randomly hit a point which is e above the stability threshold is of the order of e( N ~ 1 ' ) / 2 , so the time to determine the 
threshold with the accuracy e by the brute-force strategy grows inversely. In contrast, the method suggested herein 
allows to reach the LOCT point and determine the ST in quite a few steps. 

The estimate <[5j) can be also used to determine the probability of a wrong conclusion about the ST. Suppose we 
have a potential ST, i.e. a global minimum point Mi with ay = dist (Mi, A) and want to check if there exists another 
LOCT point M 2 with smaller ct 2 = dist (M 2 , A) = ay — 6a. For this sake we choose random starting points Ki at 
q = a± + £ until finding a point out of the basin and run the algorithm. If it converges to a point M 2 different from 
Mi and with ct 2 < a\ we have found a new potential ST. And if it converges to M\ two competing hypotheses are to 
be evaluated: 

Hypothesis A. The point Mi indeed corresponds to the ST and there is no other point M 2 . In this case, the 
algorithm converges to M\ with probability P(Mi\A) = 1. 

Hypothesis B = A. The point M 2 exists, but we have not found it. According to ©, the probabilities to converge 
to Mi and M 2 relate as 


P(M 1 \B) (, q-a 1 )^ ( 1 V V 

p(m 2 \b) a ~ (g _ CT2) ^ [l + fj 


( 6 ) 


According to Bayes’ theorem, the posterior probability P(B\Mi) is expressed through the prior probability P(B ) 
as 


P(B |Mi) 


P(Mi\B)P(B) 

P(M 1 \A)P(A) + P(M l \B)P(B) 
1 



(7) 


In Fig. 4b the posterior probability is plotted versus the prior probability for several values of a. Note that if a is 
small enough and P(B) is not close to one 0 can be approximated as P(B\Mi) « aP(B), which means that each 
trial decreases the probability by the factor a. Thus, after M trials the posterior probability of B can be estimated 
as 


In-P m 


M(N - 1) 

2 n 



( 8 ) 


The estimate ® gives the probability that after M trials we have still missed a LOCT point whose distance to the 
attractor is by Oct lower than cti. By sufficient increasing of M this mistake probability can be made arbitrarily small. 


V. CONCLUSIONS AND DISCUSSION 

To conclude, we have introduced a novel measure to describe stability of dynamical systems against external 
perturbations. This is the stability threshold (ST) which equals the magnitude of the weakest perturbation capable 
to disrupt the established dynamical regime. The ST provides important information, since it guarantees the system 
to withstand any perturbation of smaller magnitude. In the phase space, the ST is the minimal distance between the 
system’s attractor and the border of its attraction basin. From this prospective, the ST defines the “thinnest site” of 
the basin. And as the saying goes, where something is thin, that is where it tears: the direction corresponding to ST 
is the most dangerous for the system. 

For dynamical networks, different directions in the multidimensional phase space are associated with different nodes. 
To this end, the ST approach allows to determine the nodes which are mostly susceptible to perturbations. Applying 
external perturbations to these nodes, one may disrupt the network comparatively easily. However, sometimes the 
ST is associated with perturbations involving several nodes. An example of such a situation is depicted in Fig. 3a. 









Under such circumstances, it is easier to disrupt the network by simultaneous perturbation of several nodes rather 
than by perturbing just one of them. 

We have also suggested an algorithm to calculate the ST for arbitrary dynamical systems and demonstrated its 
effectiveness. The generality of the ST-based approach defines its vast potential for applications. Possible fields 
include engineering, neuroscience, power grids, Earth science and many others where robustness of complex systems 
against large perturbations is important. 

Besides application to the study of particular systems, further development and extension of the approach per 
se are of great interest. To this end, several directions are seen. The most obvious modification of the suggested 
approach is to apply it to problems opposite to the one considered herein. Particularly, in many applications the task 
is to change the behavior of the system by an external action m- From the nonlinear dynamics prospective, this 
means pushing the system from one attractor to the basin of another. In this situation, quantification of the stability 
threshold immediately provides the information about the optimal perturbation to induce the switching. 

Another possibility is to extend the concept of the stability threshold to multistable systems and reversible transi¬ 
tions between different attractors. In many applications, not a single, but several different stable regimes play essential 
roles, and the system may switch these regimes. The examples are switching between various patterns of activity in 
neural circuits mini, and episodic emergence of El Nino events m, transitions between free flow and congestion in 
traffic m or between failure states and active states in self-recovering networks [2D]. In all these cases, transitions in 
both directions are of interest. To study such transitions from one attractor to another and backwards, two stability 
thresholds can be introduced associated with each of the transitions. The values of these thresholds may provide 
important information about the transition rates when the system is externally perturbed. 

Finally, the numerical algorithm can be further developed and extended for a broader class of systems. In particular, 
the current version of the algorithm relies on smoothness of attraction basins borders, which is a limitation of the 
method. It is known that basin boundaries for some dynamical systems can be fractal |21j . In this case the border 
can not be locally approximated by a tangential hyperplane, which is crucial to determine the direction in which to 
move along it. However, it is still possible to move along the border if the steps are taken in random direction. Each 
step is accepted if it brings the point closer to the attractor and declined otherwise. This modification of the method 
requires additional investigation to the define the conditions of its applicability and estimate its convergence. 
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